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Abstract 

We study inhomogeneous host-pathogen dynamics to model the global amphibian pop¬ 
ulation extinction in a lake basin system. The lake basin system is modeled as quenched 
disorder. In this model we show that once the pathogen arrives at the lake basin it spreads 
from one lake to another, eventually spreading to the entire lake basin system in a wave like 
pattern. The extinction time has been found to depend on the steady state host population 
and pathogen growth rate. Linear estimate of the extinction time is computed. The steady 
state host population shows a threshold behavior in the interaction strength for a given 
growth rate. 


1 Introduction 

The problem of population extinction in ecological systems has been studied extensively both 
theoretically and experimentally. Population extinction can occur due to a number of reasons 
such as habitat loss, climate change, pollution, epidemics, etc m- Also, when a species alien to 
a stable ecological system is introduced it may pose a threat of extinction for a certain number 
of species [Ml- 

In an ecological system the underlying dynamics could be quite complex due to its large 
degrees of freedom m- A theoretical model is often helpful in understanding the dynamics of 
extinction without going into finer details of the real system. One approach would be to include in 
the model only the relevant degrees of freedom and incorporate the effect of the rest in parameters 
that can be determined experimentally. Population model with spatial inhomogeneity have been 
studied by a number of authors mm- Modelling extinction in inhomogeneous system can 
further our understanding of the dynamics of real systems. For instance spatial inhomogeneity 
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has a stabilizing effect on the interaction among different species jl9|B0j. it can affect the front 
speed m and the width [22] of a propagating invasion. 

The simplest example of inhomogeneous reaction-diffusion system is the trapping of diffusing 
particle by traps [23]. The density of particles for trapping reaction with randomly distributed 
traps has been know to show a stretched exponential behavior |24fl26j . self-segregation |27fl30j . 
and self-organization around traps |3H434j . Random growth and trapping of diffusing particle 
can be used to model inhomogeneous population models [ 35][36j . Recently trapping reaction 
models has been used to study population dynamics with localized predation [ 37] . 

In this paper we study host-pathogen dynamics in a spatially inhomogeneous system. Our 
motivation come from the recent studies of global amphibian population extinction due to infec¬ 
tion by the fungal pathogen Batrachochytrium dendrobatidis(Bd ) [5][6j . Detailed investigation of 
the frog-Bd dynamics in lake basins have shown that after the arrival of the pathogen in a basin 
it spreads in a wave like pattern infecting the entire population. It is likely that the spread 
of pathogen occurs via an unknown vector since the inter-lake frog movement has not been 
found. Complete extinction is observed only when the intensity of infection reached a critical 
threshold [5] . Here we study a model of host-pathogen dynamics that explains these observed 
phenomena qualitatively. We shall show how the pathogen spreads from one lake to another by 
a wave like pattern and estimate the linear wave speed. The extinction time at a lake depends 
on the steady state population of a lake. We compute the extinction time analytically and nu¬ 
merically. By numerical computations we shall investigate the host-pathogen dynamics in a lake 
basin system to understand global population decline. Further applications and improvements 
to the model are discussed. 

2 Formulation of the model 

We consider interaction between fungal pathogen ‘Bd’ and host ‘frog’ in an inhomogeneous 
media. The inhomogeneous media consists of localized regions in space that are habitats of frogs. 
These localized regions for our case models the lake basin systems described in Ref. [5], As the 
interlake movement of frogs has not been observed, we assume that the pathogen interacts with 
the host only in these localized regions. The size of lakes are assumed negligible when compared 
to the mean separation between two nearest lakes. Therefore, we assume the lakes to be point 
like regions in space, as a result host-pathogen dynamics in the lake basin system reduces to 
a host-pathogen model in presence of quenched disorder [23]. This assumption is valid as long 
as we are interested in the depletion of the total host population in the lake basin. However, 
a generalization is possible where one can replace the point objects by extended objects. Such 
a model would be interesting only when one is interested in intra-lake dynamics with lakes of 
various shapes and sizes. 

Let us denote by 6(x, t ) the density of pathogen at position x at time t and the host population 
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in the i-th lake by fi(t). Reaction-diffusion equation for the host-pathogen dynamics is given by 


dtb(x, t ) = Dd£b(x, t ) + $(x) 6 (x, t) - pb(x, t) + eb j X ’ ^ , 

1 + 6 (x, t)/bo 

fi(t) =l[l-Oifi(t)\fi(t)-ab(xi,t)fi(t), i = 1 , 2 ,..., ( 1 ) 

where ' = d/dt, D is the diffusion constant of pathogen, the function <L(x) = pJ2i <H X — x i)/*(£) 
is the spatially varying host dependent growth rate of the pathogen, p is the pathogen decay 
rate and the term e 6 /(l + b/bo) describes the growth of pathogen in the absence of host with e 
and bo constants. Note that the choice of this particular form is to ensure a bounded growth 
of the pathogen population. At low densities the pathogen growth is exponential where as 
at high density it is constant. Such growth is seen when there is a competition of resources. 
The function <h(x) models the lakes as point like object in space with the ith lake located at 
position x = X*. The boundary condition is given by limixi^^ b(x, t) = 0 and initial conditions 
b(x, 0) = d(x — xi), fi( 0) = 1 for all i = 1,2,..., N. 

The host population fi obeys a logistic equation with a growth rate 7(1 — 6ifi ) and the inter¬ 
action term —ab(xi,t)fi describes the decay of host population due to pathogen infection. Note 
that the decoupled equations (i.e. for p = a = 0 ) have homogeneous steady states b* = 0 , e/p—1 
and f* = 0, 1 /Qi. This implies that host and the pathogen can sustain their populations sep¬ 
arately provided the nonzero steady states are stable. Let t —>• Dt, p —>• p/D, p —>• p/D, e —>• 
e/D, 7 —» 7 /D, a —>• boa/D and b —> b/bo so that we have 

d t b(x, t ) = dlb{x, t ) + $(x)6(x, t) - pb(x, t ) + 1 » 

/»(*) =7([l-0ifi(t)]fi(t) -ab(x.i,t)fi(t), i = 1 , 2 ,..., (2) 

with the boundary condition lim | a .|^ >00 b(x, t) = 0 and initial condition b(x, 0 ) = S(x — xi)/bo, 
fi( 0) = 1 for alii = 1,2,... , N. 

3 Steady state and linear stability 

Let us consider the homogeneous steady states b* = 0, e/p — 1 and f* = 0, l/0j of the decoupled 
equations. For the host, the homogeneous steady state f* = 0 is unstable and f* = 1/0* is 
stable. In other words, without the influence of pathogen each lake can be considered as a stable 
ecosystems. For the pathogen, let b — b* ~ exp(*kx+cut) be as small perturbation so that Eq. ([21) 
gives 

uj = —[ k | 2 — p + e/(l + 6 *). (3) 

Clearly, the steady state b* is stable if e/(l + b*) < p. This implies that the steady state 
b* = 0 is stable if e < p. However, if e > p there exist a range of unstable modes k such that 
—yje — p < |k| < y/e — p for which which b* = 0 becomes unstable. On the other hand, the 
steady state b* = e/p — 1 is stable for e > p. 
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One would be interested to know how these steady states change when host-pathogen inter¬ 
action is introduced. We expect the host population to decline due to the pathogen interaction. 
Let us assume that a small patch of pathogen have just arrived at a particular lake i at time 
t = 0. We want to understand how this small patch grow in time. We shall estimate the time 
scale for the total extinction of the host in the lake. 

Assume that the lake i has a steady state population /,(f) = 1/0, for t < 0 and at time 
t = 0 we introduce a small perturbation b'(x, 0) = <5(x — x,)6 to the steady state b* = 0 where 
0<6«1. Linearizing Eq. ([2]) near b* = 0 and f* = 1 /0, we obtain 

d t b'(x, t ) = d 2 b\x, t ) + p6^ 1 5(x - Xj)i/(x, t) + (e - /x)i/(x, t), 

f'i(t ) = -7'/»'(*) - a0"V(x;, t), (4) 

with boundary conditions b' = 0 as |x| —>• oo and initial conditions b'(x, 0) = <5(x—x,)6, //(()) = 0. 
Note that in Eq. Q we have ignored all lakes other than the lake i as we are interested only in 
the local dynamics. We notice from Eq. Q that the equation for the pathogen can be solved 
exactly and the solution can then be used to obtain the host population /,. In order to find 
a solution b' let us transform b ' —>• b' exp(—(e — p)t) so that we have the simplified equation 
dtb'(x.,t) = d^b'(pc,t) + p9~ 1 5(x — Xj)6'(x, t). Taking Laplace transform we obtain the solution 


b'(xi,s) 


bG(xi\xj) 

1 - 9~ 1 pG(■Xi\xi) , 


( 5 ) 


where G(x|y) denotes the Green’s function defined by G = (s — c^) 1 • Substituting for the 
Green’s function, G(x|y) = exp(— v^lx — y\)/2y/s, taking inverse Laplace transform and multi¬ 
plying the factor exp((e — p)t) we obtain 




b 



+ 


_P_ 

49, 




^erfc( 



( 6 ) 


In Fig. |T| we plot the ratio b'{x.i,t)/b as a function of time. The divergence at t = 0 is due 
to the initial condition we have chosen. Initially there is a drop in the population b' due to 
the f -1 / 2 term that arises as a result of pathogen diffusing out of the lake. We observe that 
the population grows exponentially at a rate K = p 2 /49f . The fact that K is proportional to 
the square of steady state host population 1/0,; indicates that pathogen growth rates in thickly 
populated lakes will be large. 

The host population fi(t) can be written as 

fi(t) = j- - 7T / 6'(xi,r)e- 7(t - r) (iT (7) 

"z "z Jo 

Substituting the expression for h'(xj,t) in Eq. O we obtain 


1 abe 

f.(t) = — _—_ 

9i 2 9i(K + R) 


K (e 


(. K+R)t 


1 + erf(V Kt) — 1^) + Vr erfi {VRt 


( 8 ) 
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Figure 1: Growth of pathogen population in the ith lake. At t = 0 the lake is infected by 
introducing a small perturbation to the steady state. Initially it decreases as t -1 / 2 then there is 
an exponential growth. Parameters: p/Oi = 1, e — p = 1. 


with R = e — p + 7 . Now we can estimate the time scale t ext for the complete extinction of 
the host population by solving the equation fi(t) = 0. Taking only the fastest growing term in 
Eq. ([HD we obtain 


^ext 


K + R — 7 


log 


V 3 bay/K ) 


( 9 ) 


Note that to obtain Eq. Q the erf(V Kt ) has been replaced by its mean. Expression for text hi the 
original variables can be obtained by inverting the transformation we have used. The actual time 
for extinction shall however be larger than that we have obtained in Eq. Q from the linearized 
equations. The linear estimate provides only a clue as to how it will depend on the parameters. 
For example, the extinction time t ex t in Eq. ([9]) shows a sigmoidal behavior Oj. This implies that 
lakes that are sparsely populated will take a very long time for extinction. On the other hand for 
thickly populated lakes it is the opposite. In Fig. [2]we have computed numerically the extinction 
time for different values of p. We define extinction time t ext to be the time required for the host 
population to reduce to a fraction / of it initial population. Here we assume / = 0.01. We note 
that for Oi close to zero i.e. for thickly populated lake the extinction time is very small. As 
0i increases it rapidly grows for a small range of 6{ and then shows saturation. The numerical 
values obtained from Eq. ([9]) is approximately one order of magnitude less. This is expected 
since the linearized equation assumes a constant pathogen growth rate K. We know that the 
pathogen growth rate decreases as the host population decreases therefore the linear estimates 
should be regarded as the lower bound of the extinction time. The steady state host population 
f ss as a function of growth rate 7 and interaction strength a is shown in Fig. [3j For a = 0 we 
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Figure 2: Extinction time t ex t as a function of 0,. Parameters: D = 0.1, 7 = 1 , a = 0.5, // = 
0.25, e = 0.8, b = 0.1. 


know that f ss = l/0i which is independent of the growth rate 7 . When there is host-pathogen 
interaction i.e. a / 0, the steady-state changes. We choose a range of values of a and 7 keeping 
other parameters constant. We observe that for a given value of growth rate 7 there exists 
a critical value a c such that f ss vanishes for all a > a c . In the real situation each lake can 
have a different growth rate 7 j and interaction strength a* which can give rise to more complex 
phenomena. 


4 Host-pathogen dynamics in a two lake system 

To see how the pathogen spreads let us consider a two lake system in one spatial dimension 
with the lakes located at points x = 0 and x = L. At time t = 0 we assume that the lake at 
the origin is infected and the lake at x = L is free of pathogen. The pathogen spreads into the 
neighborhood of the lake by a travelling wave. To determine the speed of the travelling wave let 
us consider the one dimensional case without interaction i.e. = 0. We have 


dfb = d^b — [ib + e 6 /(l + b ). 


( 10 ) 


A travelling wave solution b(x, t ) = u(x — ct ) satisfies the Eq. (11011 so that in a moving frame we 
can write 


d?u du 
~dz +C ~dz + 


1 + u 


— (jl\u = 0 , 


( 11 ) 
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Figure 3: Steady-state host population f ss as a function of growth 7 and pathogen interaction 
strength a. Parameters: D = 0.1, p = 1.0, p = 0.25, e = 0.8, 6 = 1.0. 


with z = x — ct and boundary conditions 


lim u(z) = e/p, — 1 and lim u(z) = 0 , ( 12 ) 

z—>— OO Z—> OO 


where u = e/p — 1 and 0 are the homogeneous steady states of Eq. (fTOl) . The constant c is the 
speed of the travelling wave. Defining v = du/dz we can rewrite Eq. dill) as 


du _ dv 
dz ’ dz 


—cv — 


1 + u 


— p ) u = 0 . 


(13) 


The steady state solution of Eq. (fl3l) are (u*,v*) = (0,0) and (e/p — 1,0). Linearizing around 
the steady states we can write dX/dz = MX where X = ( u' v') T denotes a perturbation to 
(u* v*) T and the matrix 


M 


0 

e«* __ £ 1 ,, 

(i+u*y i+u* " r p 



(14) 


The eigenvalues of M characterizes the dynamics near the steady state (u*,v*). The eigenvalues 
are 


A = —c/2 ± \Jp — e + c 2 /4, for (0,0), 

A = -c /2 ± sj(e — p)p/e + c 2 /4, for (e/p — 1 , 0 ). (15) 
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Figure 4: Host-pathogen dynamics in a two lake system. Here Lake 1 and Lake 2 are the lake 
at the origin and at x = L respectively Parameters: D = 0.1, e = 0.8, p = 0.25, a = 0.5, p = 
1, 7 = l, 0 = 0.5. 

We note that (0,0) is a stable node if c > 2y/e — p = c m i n otherwise it is a stable spiral. The 
steady state (e/p, — 1,0) is a saddle node. Travelling wave moving at speed less then c m ; n are 
unphysical Em 

Now let us consider the two lake system. The lake at the origin is infected at time t = 0 
where as the lake at x = L is free of pathogen. The time required for the pathogen to reach 
the lake is t; n f < L/(2y/e — p). In Fig. Q] we have plotted the pathogen density for the two lake 
system. The lake 1 is located at the origin and lake 2 at x = 8. Initially only lake 1 is infected 
by the pathogen so we have taken a delta function initial condition centered at the origin. As 
the system evolves in time the pathogen population spreads into the neighboring region. At 
time t = 7.8 we see a small peak at lake 2 this is approximately the time when the lake is 
infected. The density of pathogen in the region between the lakes is negligible and the peak at 
lake 2 indicates that there is a sharp growth of pathogen after it has reached the lake. The host 
population therefore decreases exponentially (see Fig. [I] inset) due to the infection. Substituting 
the values D = 0.1, e = 0.8, p = 0.25, L = 8 we obtain L/c m ; n = 17.05 which is approximately 
twice of t- m f obtained numerically. 

5 Numerical result for multiple lakes 

Numerical computation for the lake basin system in two dimensions is done by finite difference 
method. The lake basin system consists of N lakes at random positions xi,... ,xtv that are 
drawn from a uniform distribution. For our case we take N = 20 and assume xi as the origin 














which is the lake infected at time t = 0. Initial host population /,(0) = 1 for all i = 1 ,N. 
We assume 0* = 0.1 for all 1 < i < 20 and p = 7 = 1. The parameters e = 0.8, p = 0.5 and 
a = 0.5 are chosen. As larger values of diffusion constant flattens the density profile 6(x, t) we 
choose as small D = 0.25 for which the travelling wave is clearly observable. In Fig. [5)^ a) we 
can see the density b(x,t) at t = 0.01 which shows a peak at the origin i.e.. lake 1. As time 
progresses the pathogen spreads out infecting the neighboring lakes. The infection spreads from 
one lake to the nearest neighboring lake. The new lakes that are infected infected acts as new 
source of pathogen which then infect their neighboring lakes. The pathogen population grows 
as long as the lakes are populated by the host. After complete extinction at a lake the pathogen 
population reaches it steady state. In Fig. HDjl) we observe that pathogen density close to the 
origin has reduced as compared to Fig. [5{b). We also observe that new lakes that are further 
away from the lakes at the origin has been infected. A wave like pattern can be seen clearly 
from Fig. [5ja-d) . We have computed the global host population as a function of time which 
is plotted in Fig. [ 6 j We observe that the global population decline qualitatively agrees with 
the experimental observations. We would like to emphasize here that quantitative prediction of 
the model can be tested only when growth and decay rates of the pathogen and the interaction 
strengths are provided. 

6 Concluding remarks 

We studied an inhomogeneous host-pathogen model to understand the global amphibian pop¬ 
ulation extinction. We assumed that the pathogen has a bounded growth and a linear decay. 
The host population on the other hand is confined to the lakes and is assumed to have logis¬ 
tic growth. Inhomogeneity is introduced as quenched disorder to model the lake basin system. 
The host-pathogen interaction is assumed bilinear in host population and the pathogen density. 
With this minimal model we have shown qualitatively the phenomena that has been observed 
experimentally. We observed that the pathogen spreads in a wave like pattern infecting one lake 
after another and eventually spreading to the entire lake basin. By linear analysis we obtained 
the extinction time and minimum wave speed which suggests that lakes that are thickly popu¬ 
lated become extinct at a higher rate. The extinction rate is proportional to the square of the 
steady state host population. We calculated a lower bound of the extinction time. 

A natural response of the system to the pathogen invasion should be the host migration. 
It would be interesting to see how extinction could be suppressed by including host migration. 
Also, to include intra lake dynamics the point like objects can be replaced by extended objects 
to model lakes of different shapes and sizes. The parameters needed in this model may be 
obtained locally from a sample of lakes assuming that all the lakes provide similar habitat for 
the host. Furthermore, for quantitative predictions parameters such as diffusion constant and 
growth(decay) rates should be determined from experiments. 
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Figure 5: Pathogen density shows wave like pattern in a multiple lake system. Here lake 1 is 
located at the origin which is infected at t = 0 and lake 2 -N with N = 20 are located at random 
positions. Parameters: D = 0.25, e = 0.8, fj, = 0.5, p = 1.0,7 = 1, a = 0.5 and 9{ = 0.1 for all 

i = l,...N. 
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Figure 6: Global population decline in a multiple lake system. Parameters: D = 0.25, e = 
0.8, n = 0.5, p = 1.0 ,7 = l,a = 0.5 and Qj = 0.1 for alH = 1,... IV. 
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